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ABSTRACT 

As a first step to a more complete understanding of the local physical processes which 
determine star formation rates (SFRs) in the interstellar medium (ISM), we have per- 
formed controlled numerical experiments consisting of hydrodynamical simulations 
of a kilo-parsec scale, periodic, highly supersonic and "turbulent" three-dimensional 

• flow. Using simple but physically motivated recipes for identifying star forming re- 
I gions, we convert gas into stars which we follow self-consistently as they impact their 

ff^ . surroundings through supernovae and stellar winds. We investigate how various pro- 

cesses (turbulence, radiative cooling, self-gravity, and supernovae feedback) structure 
the ISM, determine its energetics, and consequently affect its SFR. We find that the 
\l ■ one-point statistical measurement captured by the probability density function (PDF) 

' is sensitive to the simulated physics. The PDF is consistent with a log-normal distribu- 

fH ^ tion for the runs which remove gas for star formation and have radiative cooling, but 

^ i' implement neither supernovae feedback nor self-gravity. In this case, the dispersion, 

(Js, of the log-normal decays with time and scales with + (Mrms/2)2) where 

^ ■ Mrms is the root-mean-squared Mach number of the simulation volume, s = In p , and 

^ ' yO is the gas density. With the addition of self-gravity, the log-normal consistently 

, under-predicts the high density end of the PDF which approaches a power law. With 

supernovae feedback, regardless of whether we consider self-gravity or not, the PDF be- 
comes markedly bimodal with most of the simulation volume occupied by low density 
gas. Aside from its effect on the density structure of the medium, including self-gravity 

• and/or supernovae feedback changes the dynamics of the medium by halting the de- 
I cay of the kinetic energy. Since we find that the SFR depends most strongly on the 

underlying velocity field, the SFR declines in the runs lacking a means to sustain the 
kinetic energy, and the subsequent high density constrasts. This strong dependence 
on the gas velocity dispersion is in agreement with Silk's formula for the SFR (Silk 
2001) which also takes the hot gas porosity, and the average gas density as important 
parameters. Measuring the porosity of the hot gas for the runs with supernovae feed- 
back, we compare Silk's model for the SFR to our measured SFR and find agreement 
to better than a factor two. 

Key words: galaxies - ISM: theory ~ techniques: numerical hydrodynamics: star 
formation 



1 INTRODUCTION 

The correlation in galaxies between the star formation rate 
and the average gas surface density over several orders of 
magnitude (Kennicutt 1998) suggests a simple, determinis- 
tic prescription (Schmidt 1959) for star formation. Yet the 
finding that, at least in the Milky Way, all star formation 
occurs in dense, cold clouds of molecular hydrogen and dust 
raises the question of how information about the average 
gas density of a galaxy reaches the small scale on which star 



formation occurs. Furthermore, observations of our own in- 
terstellar medium (ISM) as well as that of other galaxies 
reveal that far from being well described by a global quan- 
tity like the average gas density, the ISM has a spectacu- 
larly complex structure on many scales. Diffuse ionized gas 
in edge-on spirals is concentrated in webs of filaments and 
shells (Rand, Kulkarni & Hester 1990, Dettmar 1992, Fer- 
guson, Wyse, & Gallagher 1996). Atomic gas detected by 
21 cm emission in our Galaxy (Heiles 1979, 1984) as well as 
in several other spirals (Irwin 1994, Rand & van der Hulst 
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1993, Lee & Irwin 1997, King & Irwin 1997) resides in "su- 
pershells" and "worms". In maps of the nearby spirals M31 
and M33 (Brinks & Bajaja 1986, Deul & den Hartog 1990), 
it is also found to be depleted in numerous 100 pc - 1 kpc 
"holes" . Attempts to quantify this elaborate ISM structure 
are confronted with questions of identification. Structures 
are interconnected, with, for example, denser regions of gas 
embedded within filaments. Hence for example, potential 
sites of star formation cannot be picked out, without intro- 
ducing a density threshold and thereby a bias to separate 
them from the underlying density field. An alternative way 
to analyze the ISM is with Fourier transform power spectra. 
Applied to HI emission maps of the Large and Small Mag- 
ellanic Clouds, power laws over ~ 2 orders of magnitude are 
found (Stanimirovic et al. 1999, Elmegreen, Kim, Stavely- 
Smith 2001), providing another insight into the structure 
of the ISM, namely that as other observations have already 
suggested, it is likely to be turbulent. 

Clues about the energy sources for the stirring of the 
ISM come from measurements of the sizes and velocities of 
shells. In some cases stellar winds and supernovae are found 
to be adequate for creating the supershells, and HI holes. 
In other cases larger quantities of energy are demanded 
and then collisions of external clouds with the galaxies are 
invoked (Tenorio-Tagle 1981). As for the diffuse ionized 
medium, although the energy available from O stars would 
be sufficient to account for its photoionization, a well-known 
problem is that photons from the O stars cannot travel far 
from their origin without being absorbed by the molecu- 
lar clouds and HI halos surrounding them. In that case the 
photons either reach larger distances by traveling through 
photoionized conduits carved out by earlier supernovae or 
as suggested by an alternative model they are additionally 
generated in turbulent mixing layers at the interfaces be- 
tween hot and cold gas. These are ubiquitous in the ISM, 
and have been invoked as an efficient means to convert the 
thermal energy generated by shear flows to ionizing radia- 
tion (Begelman & Fabian 1990, Slavin, ShuU, & Begelman 
1993). Ultimately the energy source in the latter model is 
again the supernovae which create the hot gas. Recent X- 
ray images from Chandra map out this hot, tenuous gas, 
predicted by Spitzer (1956), above and below the galactic 
plane of disk galaxies (Wang et al. 2001). Even without a 
heat source due to its long cooling time, once it is generated 
by supernovae, such gas can persist for millions of years. Cox 
& Smith (1974) reasoned that given that OB stars occur in 
associations, it is likely that a supernovae will go off inside 
the hot cavity generated by a previous supernovae, thereby 
rejuvenating it and creating an even larger cavity. In this 
way, successive supernovae can overlap creating a network 
of tunnels. Expanding at high speed within these tunnels, 
the hot gas can move above the galactic plane where it is 
either halted by insufficient speed to escape the galactic po- 
tential, or by an encounter with a large mass of cold, high 
density gas, or by efficient mixing with cooler gas which in- 
creases its density thereby accelerating its radiative energy 
losses. 

In light of this complex environment in which star for- 
mation occurs, it is even more surprising that the Schmidt 
law is so successful. It is in the context of this complexity, 
that we undertake a study of the star formation rate in a 
multiphase ISM. We restrict ourselves to a local study of 



the ISM, namely that of a ~ 1 kpc^ region. The earliest 
local study which included supernovae feedback was done 
by Rosen & Bregman (1995) in two dimensions. They con- 
sidered a segment of a galactic disk, taking into account a 
fixed external gravitational potential, but neglecting rota- 
tional effects, self-gravity, and magnetic fields. In a three- 
dimensional model which included the effects of an external 
gravitational potential, rotation, and magnetic fields, Korpi 
et al. (1999a,b) studied a supernova driven galactic dynamo. 
Meanwhile, to investigate the disk halo interaction, Avillez 
(2000) followed the evolution of a segment of a galactic disk 
with an adaptive mesh refinement code. Unlike these stud- 
ies, ours follows self-consistently and in three dimensions 
both the gas and the stars, treating the latter as a system 
of coUisionless particles subject to gravity. Rosen & Breg- 
man (1995) followed the stellar component but treated the 
stars with the same fiuid equations used for the gas thereby 
making their flow more viscous than that expected for a col- 
lisionless system of particles. Without star particles tagged 
with their ages, Rosen & Bregman (1995) decided upon a su- 
pernovae rate for their simulation, then proceeded to set off 
supernovae with a probability of occurrence correlated to the 
stellar density. Avillez (2000) approached the issue by con- 
structing an algorithm to distinguish between isolated and 
clustered supernovae. For isolated supernovae events, Avillez 
(2000) randomly determined the positions of supernovae in 
the disk plane with rates based on observed ones. To mimic 
clustered supernovae, a percentage of the supernovae sites 
were chosen to coincide with locations where there was a 
previous supernova. In the Korpi et al. (1999a,b) implemen- 
tation there was a density criteria to determine the locations 
of isolated supernovae. In both Avillez (2000) and Korpi et 
al. (1999a, b), supernovae occuring above the disk plane were 
placed in random locations with an exponential distribution 
characterized by a scale height also adopted from observa- 
tions. Given that we are interested in the impact of super- 
novae feedback on star formation, we cannot rely on these 
methods of modeling the supernovae locations. Instead we 
require that the locations, ages, and masses of the star par- 
ticles self-consistently determine the supernovae events. A 
simple calculation shows that a star with a velocity of 10 
km/s will travel ~ 100 pc (e.g. the average size of a molec- 
ular cloud) in 10 Myr. The latter corresponds to a typical 
time delay between the birth and death of a star with M 
~ 80 Mq. In a follow-up paper we explore how our results 
change when we neglect this time delay and instead allow 
the stars to explode as supernovae immediately after their 
birth (Slyz, Devriendt, Bryan, & Silk, in preparation). Ob- 
viously a local model such as ours is of limited relevance for 
quantitative comparisons to the ISM in galaxies. As later 
detailed in section|Sl the limitations of our idealized bound- 
ary conditions and the absence in our models of an external 
gravitational potential as well as of a shear flow arising from 
rotation means that there are many fundamental questions 
that we cannot address. Nevertheless we believe that for the 
purposes of studying the non-linear interplay between star 
formation and stellar feedback, our simple model yields im- 
portant insights. 

The question we address is what physical processes reg- 
ulate the rate at which gas turns into stars in a multiphase 
ISM. In section |5| we describe the numerical method we 
use as well as the ingredients of our simulation. To model 
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the large dynamic range in densities and temperatures of a 
gaseous medium compressed by self-gravity and by shocks 
maintained by supernovae and stellar winds, a robust high- 
resolution hydrodynamical scheme proves essential. To get 
a qualitative idea about the phenomena involved, section |3| 
presents the general morphological, thermodynamical, and 
dynamical features of our simulations. A more quantitative 
analysis of the gas structure and dynamics is presented in 
section |1] where we explore changes in the gas probability 
density function and energy spectra with the addition of 
more and more physics thought to be relevant for star for- 
mation. Section|S]compares the star formation rates we mea- 
sure in our simulations to simple analytic prescriptions and 
section|S|discusses the limitations of our simulations. Finally 
our main conclusions are summarized in section Q 



2 METHOD AND INGREDIENTS OF THE 
SIMULATIONS 

Traditionally the problem with numerical simulations try- 
ing to model star formation and feedback processes is that 
the radiative losses of the hot component generated by su- 
pernovae are enourmous, even though in the absence of any 
interaction of the hot gas with the cold gas the cooling time 
of the hot gas is on the order of 100 Myr. In many cases the 
culprit is numerical diffusion which mixes cold gas into the 
hot gas more than it physically should. As a result, since the 
density of the cold gas is high, mixing even a small fraction 
of it with the low density hot gas increases the density of 
the hot component sufficiently for it to cool more efficiently 
than it should. For this reason, high resolution grid codes 
are better suited for studies of the multiphase interstellar 
medium than more diffusive particle based methods which 
require carefully constructed algorithms to circumvent arti- 
ficial cooling (e.g. Marri & White 2003). 

With this in mind, we model the evolution of gas and 
stars in a three-dimensional periodic box which is 1.28 kpc 
on a side with a grid-based scheme for the gas and a particle- 
mesh method for the stars. More specifically we have incor- 
porated the BGK hydrocode (Prendergast & Xu 1993, Slyz 
& Prendergast 1999) into Bryan's ENZO code (Bryan & 
Norman 1997, 1999) which uses a Lagrangean particle-mesh 
(PM) algorithm to follow the coUisionless stars moving in 
the gravitational potential the gas and the stars themselves 
generate. Based on gas-kinetic theory, BGK computes time- 
dependent hydrodynamical fluxes from velocity moments of 
a distribution function which is a local solution to a model of 
the coUisional Boltzmann equation, namely the BGK equa- 
tion (Bhatnagar et al. 1954). The hydrodynamics code has 
been extensively tested on discontinuous non-equilibrium 
flows (see Xu 1998 for a review) and performs well both at 
flow discontinuities and strongly rarefied regions, a criterion 
which is mandatory for ISM simulations. 

Initially the gas has constant density (pgas = 1 
atom/cm^) and temperature (Tgas = 10^ Kelvin) and simi- 
lar to the initialization in MacLow et al. (1998), its velocity 
field is drawn from a gaussian random field characterized by 
a power spectrum scaling like k"''. We truncate this veloc- 
ity power spectrum so that the field only has power on large 
scales, i.e. in modes up to k = 4. The initial Vrms is ~ 50 
km/s. Contrary to MacLow et al. 1998, we do not add veloc- 
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Figure 1. Cooling curve with vertical dotted lines overplotted to 
delineate several different temperature regimes which we consider 
in section 121 



ity perturbations at each time step to drive the 'turbulence'. 
We only impose the velocity perturbations once at the be- 
ginning of the simulation. We assume radiative cooling of 
an optically thin gas which is in coUisional ionization equi- 
librium. More specifically, our cooling function, displayed in 
figure Q is an extension of the cooling curve of Sarazin & 
White (1987) down to temperatures of Tmin = 310 K to 
account for H2 cooling using the rates given in Rosen & 
Bregman (1995). The extension to lower temperatures as- 
sumes a solar metallicity, a completely ionized gas at 8000 
K and an ionization fraction that gradually drops to 10~^ 
below 8000 K. Fitting a piecewise power law to our cooling 
curve gives: 
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if r < 310K, 

if 310K <T < 20007^ 

if 2000A' <T < SOOOif 

if 8000i4: < T < 39811_K' 

if 39811A: <T < 10'" K 

if W^K <T < 2.884 X 10^ AT 

if 2.884 X 10^ K <T < 4.732 x lO^K 

if 4.732 X 10^ -ft: <T < 2.113 x lO'^K 

if 2.113 X lO^K <T < 3.981 x 10^ 

if 3.981 X W^K <T < 1.995 x lO^A" 

if T > 1.995 X m'^K 



The lower temperature cutoff of the cooling function at 
310 K is unphysical, although Rosen, Bregman, & Norman 
(1993) argue that truncating it there is a way to model the 
contribution to the ISM pressure from sources such as mag- 
netic fields and cosmic rays, which do not decrease as the 
gas radiatively cools. 



2.1 Implementation of star formation and 
feedback 

Following Cen and Ostriker (1992), we assume that star for- 
mation is inevitable if a region is contracting (V ■ w < 0), 
cooling rapidly (tcooi < tdyn and Tgas < Tmin), and is over- 
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dense (p > pcrit)- Since we check the grid on a cell by cell 
basis to see if these conditions are met, each timescale is 
computed for each grid cell. Here tdyn is the dynamical col- 
lapse timescale, i.e. fdyn ~ ■y/3.07r/(32Gptot) where ptot is 
the sum of the gas density, p, and the stellar density, tcooi 
is the cooling timescale, i.e. tcooi = kT/nA, where n is the 
gas particle number density. Tmin is the minimum of our 
cooling curve, 310 K, and pcrit for the different simulations 
is specified in tabled If all our star forming criteria are met 
within a grid cell then we convert the following amount of 
gas, Amgas = ej^^PgasAa;^ into a "star particle", where e is 
a star formation efficiency whose value is given in table ref- 
sims, and is the updating timestep. We only allow at 
maximum 90% of the gas in a cell to be converted to stars 
in one timestep. In practice however, once supernovae inject 
hot gas into the medium, the updating timestep is short as 
it is set by the hot, low density gas. As a result At < tdyn, 
and this 90% threshold is never reached. We give the new 
star particle the same velocity as the gas out of which it 
formed and we follow the stars dynamically. The star parti- 
cle is labeled with its mass, m*, its formation time, tsF, and 
the dynamical time, tdyn, of the gas out of which it formed. 

For the purposes of the feedback however, rather than 
assume that the "star particle" formed instantaneously at 
tsF, we spread the star formation over several dynamical 
times by computing the amount of gas mass that would form 
stars after time tsF to be: 

Amstars(t) = exp^ (1) 

r r 

where r = max(tdyn, 10 Myr). With this time-dependent 
star formation rate, stars form at an exponentially decreas- 
ing rate after a dynamical time. If the dynamical timescale 
of the gas in a star-forming cell is shorter than the typical 
lifespan of a massive star, i.e. 10 Myr, then 10 Myr is used in 
place of tdyn in equation^for the value of r. Then, as a crude 
model of a stellar wind, we return 25% of Amstars to the gas, 
and since this returned mass has the velocity of the "star 
particle" we alter the momentum of the gas appropriately. 
Finally assuming only the occurence of Type II supernovae, 
we add 10~^ of the rest-mass energy of Amatars to the gas' 
thermal energy (Ostriker & Cowie 1981). The supernovae in- 
put is added locally into one cell. We explore the limitations 
of our supernovae implementation in future work. As we do 
not have the resolution to follow every individual star and 
to therefore sample a realistic Initial Mass Function (IMF) 
for them, each star particle is more like a small star cluster 
with a typical mass in the range ~ 120 — 22OM0. 

Table presents the simulations we ran, listing the 
values of the parameters for star formation and feedback. 
Although we performed several simulations with a density 
threshold for star formation, pcrit, set to 1 atom/cm^ (runs 
B5, C5 and C6), for the remainder of the paper we focus 
only on the runs with pcrit = 10 atom/cm'^. This is because 
we found that dropping the density threshold by one order 
of magnitude to 1 atom/cm^ did not change the SFR by 
a factor ten, but merely by about 10% at the peak of star 
formation. As Tabled indicates, we also experimented with 
the value of t and found that taking a value of e = 0.01 (ten 
times smaller than our fiducial value) left the conclusions 
presented in this paper unchanged, i.e. the medium became 
porous and the SFR peaked at roughly the same value al- 



though with a slight time delay compared to the run with e 
= 0.1 . 



3 GENERAL FEATURES OF THE 
MULTIPHASE MEDIUM 

We begin by showing the time evolution of one of our simu- 
lations, namely B4, which includes all the physical processes 
we considered, namely "turbulent" initial conditions (as de- 
fined in section |5J , radiative cooling, self-gravity, star for- 
mation and feedback. In figure |5| we show the gas density, 
temperature and pressure in a 12.8 pc x 1.28 kpc x 1.28 
kpc slice of this run. Due to the compression caused by tur- 
bulence and self- gravity, the gas in certain regions, satisfies 
our criteria for star formation. Following their formation, 
this first generation of stars soon explodes as supernovae, 
releasing hot gas into the interstellar medium. The mor- 
phologies of the hot bubbles are extremely non-spherical due 
to the fact that the supernovae are releasing their thermal 
energy into a spatially inhomogeneous and non-stationary 
medium. Because this hot, low density gas has a long cool- 
ing time and because the star formation rate is sufficiently 
high, subsequent generations of supernovae bubbles overlap, 
filling more and more of the volume. Ultimately the density 
and temperature span more than six orders of magnitude 
in such a simulation and are anti-correlated: high density 
regions are cold, and low density regions are hot. As the 
third column in figure |21 shows, this anti-correlation results 
in near pressure equilibrium between these two phases of the 
gas. Nevertheless the dense gas is about one order of mag- 
nitude lower in pressure than the low density gas indicating 
that a thermal instability is active. Other regions which are 
out of pressure equilibrium by 1 - 2 orders of magnitude 
are those which have just experienced thermal energy input 
from supernovae. Self-gravitating gas would also appear out 
of pressure equilibrium, something we see in later stages of 
the simulation. 

The dynamical state of the stars and of the gas in dif- 
ferent temperature regimes in the simulation is summarized 
by a plot of the average velocity dispersions (fig. El • Guided 
by some of the features in the cooling curve (see figure 0, 
we divide the temperature into the following four categories: 
(I) T < 2000 K, (II) 2000 K < T < 10^^ K, (111)10^ K < T < 

4 X lO'' K, (IV) 4 X 10'' K < T. We compute the average ve- 
locity dispersion of the gas in each of these 4 regimes, and in 
addition, we calculate the mass-weighted velocity dispersion 
of the g well as the mass-weighted velocity dispersion 
of the stars. As the stars are assigned the velocity of their 
progenitor gas at formation, their velocity dispersion closely 
follows the velocity dispersion of the cold gas. Furthermore, 
we find that with the exception of the hottest phase (IV), 
the velocity dispersion of the other phases approximately 
settles to the following values: (I) 15 km/s, (II) 30 km/s 
and (III) 75 km/s. What is very striking in the plot of the 
velocity dispersions, is the high velocities (~ 500 km/s) at- 
tained by the hot, low density component of the gas. The 
densest structures which provide the raw material for star 
formation, collide and break apart, but are also subject to 
stripping via hydrodynamical and thermal instabilities when 
this hot, low density material fiows rapidly past them. The 
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Table 1. Summary of the performed runs. All of the runs have radiative cooling. The first three columns indicate whether self-gravity, 
star formation and/or feedback are activated. Pcrit is the density threshold for star formation, e is the star formation efficiency and the 
final column indicates the grid resolution. Each simulation cube is 1.28 kpc per side. 
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Figure 3. Time evolution of the logarithm of the velocity disper- 
sion in run B4 for the gas in different temperature regimes: (I) 
T < 2000 K (triangles), (II) 2000 K < T < 10^ K (plus signs), 
(III)105 K < T < 4 X 10*^ K (squares), (IV) 4 X 10*^ K < T (di- 
amonds). The crosses mark the average mass-weighted velocity 
dispersion of the gas and the asterices that of the stars. 



picture of a "violent interstellar medium" (McCray & Snow 
1979) emerges. 

Regarding the evolution of the thermal state of the gas, 
this is well portrayed in phase diagrams of the gas (bottom 
row of figure|lll which show the distribution of the mass frac- 
tion of the gas as a function of its temperature and density. 
Given our initial conditions of uniform density and temper- 
ature, if we were to plot a phase diagram of the gas at time 
t = Myr, all the gas would occupy a single point. Be- 
cause the initial temperature (10^ K) of the gas coincides 
with the peak of the cooling curve, by 9 Myr (first panel 
of bottom row of figure |3J the majority of the gas quickly 
radiatively cools to an approximately isothermal state at a 
temperature corresponding to the minimum of the cooling 



curve, i.e. 310 K. As we instantaneously imprint a spectrum 
of velocity perturbations at the beginning of the simulation, 
the gas acquires a range of density values and therefore has 
a spread in densities by this time. Thereafter, with the in- 
jection of hot gas into the medium, a tail of low density, hot 
gas appears. However as gas with temperatures 10^ K < T 
< 4 X 10® K (phase III) is thermally unstable, it gradually 
vanishes from the medium, dividing the gas into two parts 
in the phase diagram. The majority of the coldest (T ~ 300 
K) gas differs by approximately a one order of magnitude 
pressure jump from gas with T > 5 x 10^ K. Finally the 
pressure of both the hot and cold gas changes with time. It 
rises as more and more hot gas fills the simulation volume, 
a situation that would probably be difi'erent if hot gas were 
allowed to escape the box. 

Although complex, pictures of the gas density and tem- 
perature distribution in a two-dimensional slice through the 
simulation volume, do not capture the intricacy of the three- 
dimensional structure. In an attempt to display this struc- 
ture, in figure |^ we plot isodensity surfaces of the gas for 
p = 10~^, 1, 10, and 50 atoms/cm^ at 50 Myrs. It is clear 
from these figures that the hot, low density component fills 
most of the volume, while the densest regions fill the smallest 
fraction of the space, and are scattered throughout the box. 
A three-dimensional rendering of the stellar density at the 
same time instant (fig-llj, reveals traces of the imprint of the 
high density gas distribution and encouragingly bears some 
qualitative resemblance to the distribution of Ha emission 
in disk galaxies (e.g. NGC 4631, Wang et al. 2001). 



4 QUANTIFYING THE STRUCTURE AND 
ENERGETICS OF THE MULTIPHASE 
MEDIUM 

In an effort to assess what determines star formation rates, 
we systematically examine how different physical processes 
change the structure and the energetics of the interstellar 
medium. The sequence of runs listed in tableware designed 
to isolate the effects of successively more complicated phys- 
ical processes. A plot comparing the star formation rates 
for this sequence of runs (figure |7J invites us to study what 



6 Slyz A. et al. 




-4-3-3-1 01 234567 BO 123 + 5 6 



Figure 2. Time evolution of the logarithm of the gas density (first column), temperature (second column) and pressure (third column) 
in a 12.8 pc X 1.28 kpc X 1.28 kpc slice for run B4. 



keeps star formation at a minimum and alternatively what 
is necessary to drive it to high values. Resolution effects im- 
mediately manifest themselves in figure|7| The 64^ and 128'^ 
runs start from the same initial conditions. Preceeding star 
formation, feedback is non-existent, but self-gravity plays 
a larger role in the 64^ run where a grid cell of equivalent 
density to that in the 128'^ grid is 8 times more massive. 
Therefore in the 64^ run with only self-gravity (run C2), 
the SFR rises more rapidly at earlier times than for the 
comparable run performed on the 128'^ grid (run B2). Once 
feedback occurs, a mechanism supplementary to turbulence 
exists for creating high density contrasts which are stronger 



in the higher resolution runs. This causes higher peaks of 
SFR in the 128^ runs with feedback (runs B3 and B4) as 
compared to the equivalent 64^ runs (C3 and C4). On the 
other hand, feedback also creates an extra source of pressure 
to fight self-gravity which explains why the C2 run leads 
to higher SFRs at earlier times than the C3 and C4 runs. 
What remains unclear without performing a simulation at 
still higher resolution is whether the indistinguishability be- 
tween the 128'^ runs with feedback regardless of whether or 
not there is self-gravity (runs B3 and B4) are a manifes- 
tation of convergence or coincidence. However, we believe 
convergence is the more probable explanation as increasing 
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Figure 4. Time evolution of the density PDF (top row) and pliase diagrams (bottom row) for run B4 (128'^ run witli star formation, 
feedback and self-gravity). In the phase diagrams, the dotted vertical (horizontal) line marks the critical density, Pcriti (temperature, 
T^crit) for star formation. Dotted diagonal lines mark lines of constant pressure, and are labeled for the t = Myr frame: P = 10^, 10^, 
10", and 10^ ke cm ^ K. Dashed diagonal lines (labeled for the t = Myr frame) mark the Jeans length: Aj =10 pc, 34 pc, 113 pc, 380 
pc and 1.28 kpc. 
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Figure 6. Isodensity surface of the stellar density for run B4 and 
Pi, = 0.1 Mo/pc^. 



the resolution tends to increase the dominance of feedback 
processes over self-gravity. More specifically, in the case of 
the 64^ runs a rise in the SFR is driven more rapidly when 
self-gravity is included. In contrast, star formation increases 
at similar rates regardless of whether self-gravity is included 
in the 128^ runs. Therefore we do not see any reason why 



this trend should be inverted by further increasing the res- 
olution. 

Before proceeding, we calculate roughly the supernovae 
rate corresponding to the measured star formation rates in 
our simulations. In our 1.28^ kpc"^ box, typical star forma- 
tion rates are SFR ~ 0.1 — 0.8MQ/yr. Scaling these values 
to a Milky Way type galaxy, where Mmw is the mass of 
gas in the Milky Way, and Mbox is the mass of gas in our 
simulation cube, 



SFR (AfMW /Mbox) ~ 100 - SOOMo/yr. 



(2) 



For a Salpeter IMF there is approximately 1 SN/200 M©, 
implying that the typical supernovae rates in our simulation 
volume are ~ 0.5 — 4 SN/yr. Furthermore, with this scaling 
to higher mass the projected gas surface density increases 
by about 4 orders of magnitude bringing both the SFRs and 
surface densities to values representative of the starburst 
regime in the Kennicutt relation (Kennicutt 1998). 

A visual examination of a 2D snapshot of the gas 
density, temperature and pressure taken at the same time 
(t = 45 Myrs) for runs including different physics is useful 
for comparing some of the consequences of the different pro- 
cesses. Figure|H|clearly shows how self-gravity, which is a ra- 
dially directed force towards regions of locally high density, 
causes high density regions to take on a more spherical ap- 
pearance. Furthermore, all the runs without feedback have 
gas with pressure spanning over ~ 6 orders of magnitude, 
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Figure 7. Time evolution of the star formation rate for a series of runs (see table differing in their physics. The left panel displays 
the results from runs CI (diamonds), C2 (squares), C3 (triangles) and C4 (asterices). The right panel displays the results from runs 
Bl (diamonds), B2 (squares), B3 (triangles) and B4 (asterices). Symbols are the measured SFRs and the dotted and dashed lines are 
analytic models from Silk (2001). 
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Figure 8. The gas density, temperature and pressure at time t = 45 Myrs in a 12.8 pc X 1.28 kpc X 1.28 kpc slice for runs A, Bl, B2 
and B4 (see table for the specifications of each of these runs). 



and a small range in temperatures compared to the run with 
feedback. The low density gas regions in the runs without 
feedback, are cold (T ~ 300 K) and are created by adiabatic 
cooling during extreme expansion in certain regions of the 
"turbulent" medium. Another feature that appears in this 
sequence of simulations is that the dense structures in the 
run with feedback are sharper due to destruction of inter- 
mediate density material by thermal and hydrodynamical 
instabilities. 

4.1 Probability Density Function of Mass Density 

A density probability distribution function (PDF) is a sim- 
ple one-dimensional statistical measure of the structure of 
a medium. In practice for simulations performed on a grid, 
PDFs are instantaneous histograms tallying the number of 
grid cells of a certain density in the simulation. Under the 
premise that stars form in high density regions, the statisti- 
cal properties of the density field, itself nonlinearly coupled 



to the velocity field, might give clues to the process of star 
formation. Efforts to uncover how the gas density organizes 
itself in media structured by different dynamical processes 
are ongoing. Vazquez-Semadeni (1994) presented a statisti- 
cal argument to show that turbulent (random), supersonic, 
compressible fiows naturally generate hierarchical structure 
without necessitating an appeal to things like fragmentation 
in a gravitationally unstable system (Hoyle 1953). In the 
limit of very high Mach numbers these flows have a pres- 
sureless behaviour and if, in addition, self- gravity is negligi- 
ble then the hydrodynamical equations are scale-invariant. 
Consequently, whatever the density in a given region, that 
region has the same probablity of producing a relative fiuc- 
tation with respect to its normalizing density, as any other 
region in the flow. Assuming that in a random flow succes- 
sive density steps are independent, the central limit theorem 
dictates that the density distribution should be log-normal. 
And indeed, Vazquez-Semadeni's (1994) two-dimensional, 
essentially isothermal (7 — 1.0001) simulations of a weakly 
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compressible (M ~ 1), turbulent flow without self-gravity 
developed log-normal density PDFs both on the large scale 
of the simulation and in subregions within the simulation. 

Subsequently, numerical experiments of three- 
dimensional, isothermal, randomly forced, supersonic 
turbulence by Padoan, Nordlund & Jones (1997) also found 
that the gas density follows a log-normal distribution, 

ppp^ ^_(l„,_<l„p»2/2,2^ (3) 

Furthermore they observed empirically that the dispersion, 
(T, of the log-normal scales with the root-mean-squared 
Mach number, Mrms, as follows: 

a^=lu(l + (^f) (4) 
or, for the case of the linear dispersion 

0"lincar ~ ^ ■ (Oj 

These dispersion relations reflect the fact that in a medium 
with higher Afrms, the gas achieves greater density contrasts. 
Passot & Vazquez-Semadeni (1998) found the same linear 
scaling relation for the isothermal case. A formal proof for 
the lognormal PDF in the case of isothermal, supersonic tur- 
bulence was provided by Nordlund & Padoan (1999) based 
on the formalism given in Pope & Ching (1993). 

Scalo et al. (1998) and Passot & Vazquez-Semadeni 
(1998) extended this work on isothermal flows by considering 
the polytropic case. Having conducted two-dimensional sim- 
ulations including various combinations of physical processes 
(e.g. self-gravity, magnetohydrodynamics, Burgers turbu- 
lence), Scalo et al. (1998) found PDFs that were more con- 
sistent with power laws than with log-normal distributions. 
Seeking to understand this result and its discrepancy with 
previous work on isothermal flows which consistently found 
lognormal distributions, Scalo et al. (1998) performed one- 
dimensional simulations of forced, supersonic, polytropic 
turbulence and uncovered a lognormal PDF for the cases 
where either the gas was isothermal (7 = 1) or where the 
Mach number was small (M^ 1). Otherwise, when 7 < 1, 
power laws developed for densities larger than the mean. 
Alternatively, Nordlund & Padoan (1999) interpreted Scalo 
et al.'s results for the PDFs occuring in the 7 5^ 1 case as 
skewed log-normals and Passot & Vazquez-Semadeni (1998) 
provided a mathematical framework for understanding why 
these distributions arose. 

Our work extends these investigations on the PDF in 
the direction of the cases where the ISM is constrained nei- 
ther to be isothermal nor polytropic. As a result our local 
temperature and pressure are not simple functions of the 
density but arise from the evolution of the thermal energy. 
Because we consider processes (e.g. radiative cooling, self- 
gravity, star formation) whose effectiveness depends on the 
density, the hydrodynamic equations are no longer scale- 
invariant. Therefore the condition of randomness between 
subsequent density fluctuations is violated and one cannot 
expect a log-normal density PDF (e.g. Vazquez-Semadeni's 
(1994)). In our series of experiments of increasing complexity 
(see Table 0, the simplest simulation we performed was of 
non-isothermal supersonic turbulence (run A). Despite the 
inclusion of density-dependent cooling processes, we found 
that the structure of the gas quickly evolved to a density 



PDF consistent with a log-normal. This is not a surprising 
result since without a heat source the majority of the gas 
quickly cools to a nearly isothermal state (see bottom row 
of flgure 1^ with an average temperature corresponding to 
the minimum of the cooling curve (horizontal dashed line 
in bottom row of flgure |5J. Furthermore, the scaling for the 
dispersion of the PDF given by Padoan, Nordlund & Jones 
(1997) continued to hold. In fact, rather than flt log-normal 
functions to our density PDFs, we measured the average 
of the logarithm of the gas density, < logj^gp >, and the 
Afrnis of the gas at different time instances and then over- 
plotted Padoan, Nordlund & Jones' (1997) prediction for 
the log-normal distribution. For the runs where we formed 
stars (without self-gravity or feedback) in addition to hav- 
ing radiative cooling (runs Bl and CI), the gas density PDF 
continued to have the same behavior: the Afrms of the sys- 
tem progressively declined with time, while the density PDF 
remained consistent with a log-normal distribution (flg. llUH . 

The runs which showed the first departure from log- 
normal density PDFs were the runs which included self- 
gravity (runs B2 and C2) but still no feedback (flgure ITTt . 
Repeating the exercise of measuring the average of the log- 
arithm of the gas density, < logjQ/9 >, and the M^-ms of 
the gas at different times, we found two differences: (a) the 
Mr-ms initially declined but then stabilized at a value higher 
than that seen in the runs without self-gravity, and (b) the 
log-normal PDF predicted by Padoan, Nordlund & Jones 
(1997) consistently underpredicted the distribution at high 
gas density. A power-law flt the high density tail well. In 
one-dimensional simulations of Burgers flows, i.e. infinitely 
compressible flows, power-laws were also found to be good 
flts to the density PDFs (Gotoh & Kraichnan 1993). We 
therefore interpret the power-law behavior for the run with 
self-gravity, as reflecting the added possibility of the gas, 
once it has a high density, to compress to even higher den- 
sity, reminiscent of the behavior in Burgers flows. Klessen 
(2000) also explored the form of the density PDF for the 
cases of decaying and driven self-gravitating turbulence. Al- 
though he found a departure from log-normal at high den- 
sities, the departure could not be characterized by a power 
law. 

When we add feedback to the list of simulated processes, 
either with self-gravity (runs B4 and C4) or without (runs 
B3 and C3), the density PDF becomes markedly bimodal 
(flgure [T^ ■ illustrating that the majority of the simulation 
volume is occupied by low density gas. A bimodal density 
distribution is also a sign of a thermal instability (Vazquez- 
Semadeni, Gazol & Scalo (2000)) the consequences of which 
we will discuss in a future paper (Slyz, Devriendt, Bryan & 
Silk, m preparation). For the runs with self-gravity, the high 
density power-law tail disappears. Perhaps it can be argued 
that the high density part of the density PDF may be flt 
with a log-normal distribution (figure [T^ . The exercise of 
overplotting the log-normal given by Padoan, Nordlund & 
Jones (1997) is not possible because the Afrms measured for 
the entire simulation box does not correspond to the Afrms of 
the high density gas for which the log-normal function may 
be a good description. Hence we can only fit log-normals to 
the high density gas, similar to what others, e.g. Wada & 
Norman (2001), Kravstov (2003), do in their global simula- 
tions of the ISM. 

The interest of describing the density structure of the 
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Figure 9. Time evolution of the density PDF (top row) and phase diagrams (bottom row) for run A (128"^ run with no star formation, 
no feedback and no self-gravity). The thick dashed line overplotted on the measured PDFs (symbols) is the log— normal PDF predicted 
by Padoan, Nordlund & Jones (1997). In the phase diagrams, the dotted vertical (horizontal) line marks the critical density, Pcriti 
(temperature, T^rit) for star formation. Dotted diagonal lines mark lines of constant pressure, and are labeled for the t = Myr frame: 
P = 10'^, 10^, 10"', and 10^ kg cm'^ K. 
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Figure 12. Comparison of the PDFs at 110 Myrs for runs with 
different physics. 



ISM with a single function, such as the log-normal, lies in 
finding a link between the gas density averaged over kilopar- 
sec sized regions and the high density regions which might 
form stars. This is precisely the link required for an explana- 
tion of the Schmidt law. Rewriting the Schmidt law in a form 
where the star formation rate is equal to some constants 
multiplied by the fraction of gas in high density regions and 
by the gas density averaged over large scales (his equation 
7), Elmegreen (2002) emphasized that star formation rates 
depend on the geometry of the density field, i.e. the PDF. 
If the shape of the density PDF is universal, then the frac- 
tion of gas in high density regions is known. Consequently, 
if the high density regions are also self-gravitating, then the 
fraction of gas available for star formation is also known. 
Admittedly, the density PDF contains no spatial informa- 
tion, hence there is no reason for which the high density 
regions should find themselves to be spatially contiguous, so 
that they comprise regions of mass greater than the Jeans 
mass. In fact, figure [TTI clearlv shows that at least some of 
the dense gas regions are not contiguous because if they were 
they would simply not persist as all the gas would be con- 
verted to stars on a dynamical timescale since these regions 
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Figure 10. Time evolution of the density PDF (top row) and phase diagrams (bottom row) for run Bl (128^ run with star formation, 
no self-gravity, and no feedback). The thick dashed line overplotted on the measured PDFs (symbols) is the log-normal PDF predicted 
by Padoan, Nordlund & Jones (1997). In the phase diagrams, the dotted vertical (horizontal) line marks the critical density, Pcriti 
(temperature, T^rit) for star formation. Dotted diagonal lines mark lines of constant pressure, and are labeled for the t = Myr frame: 
P = lO'', 10^5, 10*, and 10^ kg cm'^ K 



are well above pcrit and cold. We therefore have to identify 
these regions with divergent gas flows. 

A two-dimensional study of the ISM in a galactic disk 
by Wada & Norman (2001) has claimed that the log-normal 
distribution is a robust description of the ISM density distri- 
bution over many orders of magnitude in density, regardless 
of the simulated physics. More specifically, in their simula- 
tions the presence of stellar feedback does not change the 
shape of the PDF but increases the dispersion of the log- 
normal. In three-dimensional simulations of a high-redshift 
galaxy performed in a cosmological context, Kravtsov (2003) 
finds a density distribution similar to Wada & Norman's 
(2001). Its shape at every redshift epoch has a flat region 
at Pgas < 1 - 10 Mq pc~'^ and a power law distribution 
at high densities. He claims that the log-normal distribu- 
tion is a fair description of the high density tail of the PDF 
and agrees with Wada & Norman (2001) on the insensitiv- 
ity of the distribution to feedback, except at the low den- 
sity end, where the simulation with feedback produces more 
low density gas. As figure WH shows, our less realistic study 
of star formation occuring in a periodic box without the 
global gravitational galactic potential or the shear instabil- 
ities present in a self-gravitating rotating disk, appears to 



be more sensitive to the input physics. Only the runs which 
include stellar feedback are nearly equivalent, regardless of 
whether there is self-gravity. When log-normals are over- 
plotted for the runs without feedback, the position of the 
maximum of the log-normal is shifted to lower densities by 
more than one order of magnitude from the position of the 
maximum of the log-normal fit to the high density part of 
the PDF for the runs with feedback. Indeed the densities in 
certain cells for the run with only self-gravity reach the same 
high values as the runs with feedback, but a much smaller 
fraction of the simulation volume has these high densities. 
Another blatant difference between the PDFs we find in our 
runs with feedback and the PDFs found by Wada & Nor- 
man (2001) and Kravtsov (2003) is that their runs do not 
show as high a peak at low densities. The smaller quantity 
of low density gas in their simulations is likely due to the 
much lower supernovae rates in Wada & Norman's (2001) 
simulations (0.01 SN/yr as compared to 0.5 - 4 SN/yr in 
our simulations) and in Kravtsov's (2003) case, to the more 
realistic boundary conditions, which allow tenuous, hot gas 
to escape the disk. 
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Figure 11. Time evolution of the density PDF (top row) and phase diagrams (bottom row) for run B2 (128^ run with star formation 
and self-gravity but no feedback). The thick dashed line overplotted on the measured PDFs (symbols) is the log-normal PDF predicted 
by Padoan, Nordlund & Jones (1997). The solid line with a slope of -1.5 plotted at t = 85 Myr is a fit to the high density end of the 
PDF. In the phase diagrams, the dotted vertical (horizontal) line marks the critical density, pcriti (temperature, T^rit) for star formation. 
Dotted diagonal lines mark lines of constant pressure, and are labeled for the t = Myr frame; P = 10®, 10^, 10"*, and 10^ kg cm K 



4.2 Energy Spectra 

Energy spectra of the ISM carry complementary information 
to that given by a study of its density structure. With the 
density PDFs, we confirmed that in many cases there exists a 
clear relationship between the density contrast achieved and 
the Mrms of a system (i.e. aiinoar ~ Mrms). But the Mrms of 
a system is only a global measurement of its kinetic energy 
content. With measurements of the kinetic energy spectra, 
we expect to learn how the energy is distributed on different 
spatial scales and how the different physical processes we 
considered influence the time evolution of this distribution. 

The Kolmogorov theory of incompressible, subsonic tur- 
bulence predicts that energy fed on large scales progressively 

cascades to smaller scales until it is dissipated by molecular 
viscosity on the smallest scales in vortex rings. The transfer 
of energy is a local process and the spectra of the velocity 
field is a power law with Ek ~ k~^^^ (Kolmogorov 1941). 
With supersonic, compressible turbulence, strong shocks 
come into play. They allow energy to be transferred over 
widely separated scales and it is possible that rather than 
being dissipated in vortex rings, the energy is ultimately 
dissipated in sheets, filaments and cores (Boldyrev 2002). 



Given the analogy between highly supersonic and pressure- 
less flows, one might expect the compressible, supersonic 
flows to have the same behavior as Burgers turbulence with 
power spectra in the inertial regime of the form, Ek ~ fc~^ 
(Burgers 1974, Gotoh & Kraichnan 1993). However this ap- 
pears to only be true in one and two dimensions. In three 
dimensions, compressible, supersonic flows differ from Burg- 
ers flows because they generate vorticity (Boldyrev 2002). In 
three-dimensional simulations of compressible, supersonic, 
magnetized forced turbulence with Mach number initially ~ 
10, Boldyrev, Nordlund & Padoan (2002) find energy power 
spectra in the inertial range to be E^ ~ A;"^ '^*, i.e. close to 
the Kolmogorov value. 

As we lack the grid resolution to ascertain if the energy 
spectra in our simulations are tending towards power laws, 
we cannot make any credible statements about the values 
of the power law slopes. Furthermore in incompressible tur- 
bulence, the energy spectrum is a power law in the inertial 
regime (at k wavenumbers below the energy injection scale 
but above the energy dissipation scale). In our simulations 
the feedback energy is injected on scales equivalent to the 
grid resolution, i.e. the smallest scales, but it can propagate 
to larger scales depending on the ISM dynamics. Therefore 
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Figure 14. Time evolution of the compressible (Ec) and solenoidal (Es) components of the energy spectra for runs Bl, B2, B3 and 
B4. Symbols denote energy spectra at time intervals separated by 30 Myrs. Solid line represents time 4 = Myr. Plus signs: 30 Myr, 
asterices: 60 Myr, filled diamonds: 90 Myr, open diamonds: 120 Myr, open triangles: 150 Myr, crosses: 180 Myr, open squares: 210 Myr. 
We also draw a solid lino through the symbols when they represent the final timestep that we are displaying. Thick dashed lines indicate 
power laws with slopes similar to that of the last curve shown. 



for the runs with feedback the inertial regime has a more 
complicated meaning. Instead in figure ITU we focus on the 
time development of the energy spectra, and the presence 
of characteristic features. The standard approach involves 
dividing the kinetic energy into two components: a com- 
pressible one for which V x «comp = 0, and a solenoidal one 
with V • Usoi ~ 0. In words, the compressible component 
measures the strength of the shocks in the system, while 
the solenoidal component measures the degree of rotation. 
Typically, the compressible component is expected to decay 
faster than the solenoidal component as the shock energy is 
transformed into vortical eddy motions. 

Because we remove gas from the system to form stars, 
the kinetic energy whose spectra we measure, is rather a 
specific kinetic energy, i.e. we divide the instantaneous total 
kinetic energy by the total gas mass present at that mo- 
ment. In all our runs, the kinetic energy which is initially 
imprinted only on large scales quickly (within ~ 30 Myr) 
redistributes itself to smaller scales as well. Following this 



redistribution, for the run with neither self-gravity nor feed- 
back (run Bl), the compressible and solenoidal components 
of the energy spectra progressively decay all the while main- 
taining approximately the same form. The ratio Ec/Es is 
always less than 1, i.e. the compressible component decays 
faster than the solenoidal one, but increases towards the dis- 
sipative regime. In high resolution simulations (512'', 1024"^) 
of decaying compressible turbulence with Mach number ini- 
tially on the order of 1 (an order of magnitude lower than the 
initial Mach number in our simulations), Porter, Woodward 
& Pouquet (1998) find a similar result with Ec/Es ~ 0.1. 
In contrast to these runs in which the kinetic energy de- 
cays, the runs with self-gravity (run B2) and/or feedback 
(runs B3, B4), show energy spectra which climb to higher 
amplitudes with time and have shallower slopes than the 
decaying run (Bl). Furthermore in plots of the ratios of the 
compressible to solenoidal components, between 90 and 150 
Myrs the runs with feedback show a peak at ~ 65 pc con- 
sistent with what one would predict for the characteristic 
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Figure 13. Log-normal fit to high density end of the PDF for run 
B4 at t = 85Myr. The scaling we use for the fit is a log-normal 
with average density of 10^'^ atoms/cm^ and with a dispersion 
of lO""^'^^ atoms/cm'^. 

lengthscale for a simulation with supernovae expanding into 
a medium with ambient pressure of P = 10^ cm~^ K fcs. 
More explicitly, ignoring adiabatic and radiative losses, a 
supernovae with 10^^ ergs of energy will be halted by an 
ambient medium at this pressure when it has expanded to a 
radius, r ~ [E/PY^^ ~ 65 pc. This signature in Ec/Es for 
the run with feedback points to a way to understand SFRs, 
which we explore below. 



5 NUMERICAL VERSUS ANALYTICAL STAR 
FORMATION RATES 

An alternative to searching for a generic density PDF as 
an explanation for star formation rates, is to consider ar- 
guments concerning the competition between the expansion 
of supernovae remnants and the pressure which halts them. 
In this vain. Silk (1997, 2001) developed porosity models 
of a regulated ISM. Introduced by Cox and Smith (1974), 
porosity, Q, is proportional to the product of the super- 
novae rate per unit volume and the maximum extent of the 
4-volume of the supernovae remnants. In other terms, the 
porosity measures the fraction of hot gas, f^, in the ISM 
through the relation Q = — ln(l — fj,). Silk reasoned that 
since the supernovae production rate is proportional to the 
star formation rate (SFR), and the maximum extent of a 
supernovae remnant is limited by the ambient pressure, the 
following expression arises: 

Q = SFR G~i/' p-'/' (agas/af)-^-'^ (6) 

where pgas is the gas density, o-gas is the gas velocity dis- 
persion, and (Tf is a fiducial velocity dispersion that is pro- 
portional to i^sil^ rrigf^ (^"''■^. Here Esn is the energy of a 
single supernova, ^g is the metallicity relative to solar of the 
ambient gas, and msN is the mean mass in newly formed 
stars required to produce a supernovae. For _Esn ~ 10^^ erg, 



^g = 1, and msN = 25OM0 i.e. the case where we assume 
only the occurrence of type II supernovae with a Miller- 
Scalo IMF, the fiducial velocity dispersion is ~ 22 km s~^. 

Our simulations with feedback provided a laboratory 
to test this analytic description of the SFR. For the purpose 
of computing the porosity of the medium, we measured the 
fraction of hot gas in our volume, defining hot to be gas with 
temperature T > 4 xlO® K. For pgas in equation|S|we took 
the average gas density in our simulation volume, and for 
(Tgas we took the average mass-weighted velocity dispersion 
of the gas. We kept the value for erf at 22 km s~^ . Given these 
values as functions of time, we plotted as dotted and dashed 
lines the expectation from eq. |S| for the SFRs in figure |7| 
Computing the actual star formation rates in the box by 
defining the mass of newly formed stars to be the mass of 
stars formed in the past 3 Myr, we overplotted the results 
as symbols in the same figure. Astonishingly, the analytic 
values match the measured rates to better than a factor 2. 

Given the simplifications in the derivation of the ana- 
lytic model, there was no o priori reason for the fit to be 
a good description of the star formation rate in an inho- 
mogenous, non-stationary model of the ISM. For example. 
Silk takes the expression for the 4-volume of the SNR rem- 
nant in its cooling phase from Cioffi, Mckee & Bertschinger 
(1988). They derive it under the assumptions that the super- 
novae expands in a spherical manner, the ISM is homoge- 
nous and uniform (i.e. no density gradients), there is no 
dust cooling or thermal conduction, and the ambient ISM 
pressure is negligible until the last stage of supernovae evo- 
lution when the remnant merges with the ambient ISM. In 
contrast, we find that at least in the initial stages of our sim- 
ulations, the supernovae remnants are highly non-spherical, 
the ISM is inhomogeneous with ubiquitous density gradients 
and the ambient ISM gas pressure is highly non-negligible 
(P = 10^ — 10® cm~'^ K fcs). However, as more of the gas 
turns into stars, and the hot phase fills the majority of the 
simulation volume, the ISM does start to resemble some- 
thing more in line with the Cioffi et al. assumptions. 

When we examine in figure lT^ the time evolution of each 
of the physical quantities entering into the analytic model 
for the SFR, we find the following. The runs (Bl and CI) 
which produced the lowest star formation rates have zero 
porosity and high fractions of cold gas (/coW ~ 0.8-0.9), 
but a continuously declining velocity dispersion. The runs 
reaching a peak (runs B2, B3, B4, C2, C3) or multiple peaks 
of high star formation (run C4) all displayed depleted cold 
gas fractions after their final star formation peak, a rise to a 
maximum in its velocity dispersion at the peak, and either 
zero porosity for the case of the runs with self-gravity but 
no feedback (runs B2 and C2) or a porosity that levels off 
to a constant value around the time of the SFR peak (Q 
~ 4-5 for the 64^ case (runs C3 and C4), and Q ~ 4 for 
the 128^ case (runs B3 and B4) after the SFR peak). We 
interpret the behavior in these parameters as refiecting the 
importance of a high velocity dispersion for generating high 
SFRs. Indeed in the analytic model for the SFR (eq. |SJ|, 
the gas velocity dispersion, o-gas, plays the most important 
role, as it is raised to the highest power in the expression. 
However even with velocity dispersions sustained at high 
values ((Tgas ~ 20 km/s), SFRs will drop if the reserves of 
cold gas decline. 
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Figure 15. Plots comparing the time evolution of the cold mass 
fraction, porosity, average gas density, and mass- weighted velocity 
dispersion for runs containing different physics, on the 64^ and 
128^ grid. 



6 DISCUSSION 

Given the simplicity of our simulations, we examine their 
relevance for representing true star formation processes in 
real galaxies. The first issue we address is whether the star 
formation rates we obtain are consistent with the Kennicutt 
relation. In section |31 we scaled the mass in our simulation 
volume to that of the Milky Way, finding that our star for- 
mation rates and surface densities were consistent with star 
formation occuring in the starburst regime. If we do not scale 
our SFRs and gas densities to a Milky Way type galaxy but 
instead take them at face value we find that our initial 1 
atom/cm"^ gas density in a (1.28 kpc)"^ volume yields in pro- 
jection about a 30 M0 /pc^ column density which lies at the 
boundary between Kennicutt's normal disks and centers of 
normal disks (Kennicutt 1998). Transforming our average 
star formation rate of 0.2-0.3 Mq /yr into a star formation 
rate per unit volume leads us to an average star formation 
rate density of about 0.1 Mo/yr/kpc'^, on the high side but 
in fair agreement with Kennicutt's measurements for our 
computed surface density (Kennicutt 1998, figure 6). We 
note that Kennicutt's law is a static relation as it concerns 
space averaged quantities in local galaxies, and a moment 
in the history of these galaxies is bound to exist when their 



main progenitor will be entirely gaseous (i.e. with no stars 
yet formed) and the Kennicutt relation will break. As our 
simulations start from an exclusively gaseous medium, we 
do not expect our simulation to follow the Kennicutt rela- 
tion from the very beginning, but to move towards it as it 
does. We nevertheless consider our simulations to be in a 
starburst mode because the duration of the star formation 
episode is much shorter than that of what one expects in ei- 
ther a disk or spheroidal galaxy. But this is not unusual since 
we are only modeling a chunk of a galaxy and are therefore 
neglecting effects on larger length and therefore timescales. 

The second issue we address is whether periodic bound- 
ary conditions drive the high star formation rates seen in our 
simulations. When hot gas starts to fill the bulk of the sim- 
ulation volume, because the boundary conditions trap the 
hot gas, conditions in the simulation may be viewed as a 
pressure cooker and the increased pressure may drive higher 
star formation rates. In our simulation by the time the pres- 
sure cooker is operative, the SFRs are already at starburst 
levels as seen when one scales the SFRs and gas densities to 
a Milky Way type galaxy as we do in section |1] To be more 
specific, for the pressure cooker to be operative we have to 
wait ~ 10 Myr for the first supernovae to go off and then 
we have to wait for the volume to become significantly filled 
by this supernovae generated hot gas for the hot gas to be 
able to traverse the volume unobstructed by cold, dense gas. 
According to figure 15, it takes on the order of 50 Myr for 
the hot gas filling fraction to be approximately 50%, corre- 
sponding to a porosity of about 0.7. Hence boundary effects 
are not dominant in shaping the star formation rate un- 
til after that time. We also point out that the limitations of 
the boundary conditions should not obfuscate the point that 
the manner in which we implement supernovae is a more im- 
portant factor leading to the build up of large quantities of 
hot gas in the medium. When we perform simulations in all 
points identical to those presented in this paper but with su- 
pernovae going off instantaneously, as opposed to exploding 
with a more realistic 10 Myr time delay used in the work 
presented in this paper, we get extremely low star forma- 
tion rates (a few hundred times smaller than those we get 
in our simulation here), because the hot gas never fills a sig- 
nificant fraction of the simulation cube. In other words, the 
periodic boundary conditions cannot dominate the physics 
of star formation driven by hot gas pressure until the hot 
gas has already been generated, and we find that this de- 
pends strongly on the way the supernovae are implemented. 
As mentioned in section Q we leave the discussion of this to 
a future paper. 

The limitations of our closed, periodic box, and the ab- 
sence of a stratified external gravitational potential certainly 
keep our simulations far from being representative of realis- 
tic galactic systems. For example, a credible simulation of a 
disk galaxy, would have to be performed in a realistic cosmo- 
logical context to capture such effects as tidal encounters and 
stripping from neighbours. Excluding these external stellar 
heating processes as well as spiral waves, results in the ne- 
glect of processes that would increase the velocity dispersion 
of the stars in real galaxies. Therefore our simulations cer- 
tainly have a higher fraction of cold ISM and cold stars after 
a gas consumption time which may prolong and strengthen 
star formation in our simulations. 

We also emphasize that with our crude assumption of a 
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closed box not only can no material escape the box, affecting 
star formation rates once hot gas permeates the simulation 
volume, but no material can enter the simulation volume ei- 
ther. It could well be that accretion of cold material is more 
relevant for star formation in real disks than cither the ex- 
ternal star heating processes missing from the simulations 
discussed above or the fact that hot gas cannot leave the 
simulation volume. One can argue that perhaps the Simula^ 
tions presented in this paper are more representative of what 
happens in the central kiloparscc of a spheroidal starburst 
galaxy. In that case the potential well might indeed trap a 
fraction of the hot gas and the pressure cooker environment 
which comes into play after high star formation rates oc- 
cur in the simulation, if not as drastic as in our simulations 
might well be fairly realistic. 



7 SUMMARY AND CONCLUSIONS 

To unravel which global paxameters control star formation, 
we have examined star formation occurring in media whose 
dynamics arc structured by various combinations of physical 
processes (e.g. "turbulence", radiative cooling, self-gravity, 
feedback from supernovae and stellar winds). We sought to 
understand our models of the ISM from structural and dy- 
namical perspectives, finding that in some cases there was 
a well-defined link between the two. In particular, measure- 
ments of the density PDFs confirmed that for the simula- 
tions without feedback, lognormals were an adequate de- 
scription of the structure of the medium, and that the den- 
sity contrasts achieved in the media were directly correlated 
to their Mmis- Lognormals consistently underpredicted the 
high density end of the runs with self-gravity which appeared 
to be well-fit by a power law. For the runs with feedback, 
the dense gas reached higher densities than those reached 
by the runs without feedback implying that in these simula- 
tions, feedback was positive in the sense that it encouraged 
higher star formation rates. However the PDF for the runs 
with feedback had a distinctly bimodal shape with the ma- 
jority of the volume filled by low density gas. In summary, 
we did not find a universal PDF. Most markedly, runs with 
feedback had a different PDF from the runs without feed- 
back, although arguably, the high density end might be fit 
by a lognormal. 

Measurements of the energy spectra in our simulations 
were consistent with the information provided by the density 
PDFs. Self-gravity alone was sufficient to sustain the kinetic 
energy of the medium, and hence maintain the high density 
contrast we observed in the PDFs. Feedback also succeeded 
in keeping high quantities of kinetic energy in the media 
and inspection of ratios of compressible to solenoidal energy 
revealed that supernovae were pumping energy into the sys- 
tem at a characteristic scale consistent with the ambient 
pressure in the hot, low density component of the medium. 

For the runs with feedback, comparing Silk's (2001) star 
formation model to the measured values of the SFRs in our 
simulations, revealed a good match that led us to inspect 
the parameters involved in Silk's prescription. They showed 
clearly that the SFR depends strongly on the underlying 
velocity field which we saw could be energized by self-gravity 
and/or feedback to produce high density contrasts. Without 



a means to create these high densities, star formation rates 
decline even in the presence of a large reservoir of cold gas. 

In light of the issues neglected in our simulations, we 
stress that the simplifying assumptions made in this paper 
facilitated our choice to start from as strong as possible a 
local physical basis as possible before trying to tackle star 
formation in a more global context. As such we neglect nu- 
merous physical processes which may invalidate partially 
or completely our current results, but this remains to be 
addressed in future work. Nevertheless we hope that the 
present work sheds some light on the local physics that 
should be included in future realistic simulations of star for- 
mation. 
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